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Abstract. We present structure calculations of neutral and singly ionized Mg clusters of up to 30 atoms, 
as well as Na clusters of up to 10 atoms. The calculations have been performed using density functional 
theory (DFT) within the local (spin-) density approximation, ion cores are described by pseudopotentials. 
We have utilized a new algorithm for solving the Kohn-Sham equations that is formulated entirely in 
coordinate space and, thus, permits straightforward control of the spatial resolution. Our numerical method 
is particularly suitable for modern parallel computer architectures; we have thus been able to combine an 
unrestricted simulated annealing procedure with electronic structure calculations of high spatial resolution, 
corresponding to a plane-wave cutoff of 954eV for Mg. We report the geometric structures of the resulting 
ground-state configurations and a few low-lying isomers. The energetics and HOMO-LUMO gaps of the 
ground-state configurations are carefully examined and related to their stability properties. No evidence 
for a non-metal to metal transition in neutral and positively charged Mg clusters is found in the regime of 
ion numbers examined here. 
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1 Introduction 

The properties of nanometer-sized clusters of atoms are 
significantly different from those of the isolated chemical 
species as well as the bulk material, and generally exhibit 
a strongly non- monotonous size-dependent behavior pQ. 
Clusters are the ultimate nanostructures, where a funda- 
mental understanding of their properties can be achieved 
one atom and one electron at a time. They thus open a 
unique window to study the emergence of the properties 
of macroscopic matter from its microscopic constituents, 
and have become an active field in basic research. Nanos- 
tructured materials also hold, due to low dimensionality 
and unique composition, the promise of technological ap- 
plications reaching as far as clean and sustainable energy, 
reactions, and catalysis [2,3 . 

The transition from microscopic to macroscopic be- 
havior is particularly obvious for clusters of atoms that 

a This work was supported by the Austrian Science Fund 
FWF under projects P18134 and P21924 (to EK) and J2936 
(to SJ). Computational support was provided by the Central 
Computing Services at the Johannes Kepler Universitat Linz, 
we would especially like to thank Johann Messner for help and 
advice in using the facility. 



are bound by covalent bonds at small atom numbers, but 
are metals in bulk quantities. These must, at some point, 
undergo a transition to metallic behavior. 

The history of the examination of these clusters is long; 
we take a fresh view for a number of reasons: One is that 
we present here the first large-scale application of limerec 
[4], a new open source Density Functional Theory (DFT) 
package specifically designed for cluster calculations. The 
method is formulated entirely in real space. It thus avoids 
any basis set bias and is particularly suited for modern 
parallel computing environments [5 ,6 ,7 . Due to its effi- 
ciency, the method allows to find ground state configu- 
rations and low-lying isomers of fairly large clusters by 
unrestricted minimization of the total energy in the whole 
configuration space, using Langevin- Monte- Carlo and/or 
steepest descent annealing procedures. It thus also elimi- 
nates a possible bias due to the choice of a certain sym- 
metry or starting configuration that is present in more 
restricted minimization schemes. Our second objective is 
to prepare for studying the influence of a quantum fluid 
matrix on the formation and structure of metallic clusters 
in the near future. 

The main thrust of this paper is to present a detailed 
and systematic characterization of neutral and singly ion- 
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ized magnesium clusters of up to 30 atoms. We also dis- 
cuss some sample calculations for smaller sodium clusters, 
mostly for the purpose of studying electron localization. 
All calculations have been performed in the framework of 
spin-density functional theory (SDFT), using the Perdew- 
Wang exchange-correlation functional [8] and local [9|fTU] 
as well as norm-conserving non-local pseudopotentials [11] 
to treat the core electrons. In addition to the spatial struc- 
tures of the ground state configurations and low-lying iso- 
mers, we also report binding energies, ionization energies 
and fragmentation energies as well as the HOMO-LUMO 
gap for all clusters. Details of the computational method 
are given in section [2] 

We have chosen to study mostly Mg n and Mg+ clusters 
for a number of reasons: The simpler analog, Na n clusters, 
has been studied extensively, see, for example, Ref. [12] 
for a recent review as well as references to earlier work. 
Although there are still some open issues [13] , Na n clusters 
are reasonably well described by a simple jellium model, in 
particular energetic quantities like magic numbers are well 
reproduced. In Mg n clusters, electrons are more strongly 
localized and, hence, a jellium model is less appropriate. 
Comparing Mg n and Na n clusters can therefore tell us 
about the consequences of electron localization. 

A striking feature is that Mg n clusters evolve, with 
increasing ion number n, from molecule-like complexes 
bound by covalent bonds to a bulk metal where the valence 
electrons are delocalized. For Mg clusters in helium, this 
non-metal to metal (NMM) transition has been reported 
to occur around n = 20 with different experimental tech- 
niques [14,15,16]. Thomas et al. [16], for example, have 
measured photoelectron spectra of mass-selected magne- 
sium cluster anions, which are related to the HOMO- 
(HOMO-1) gap. Similar results have been reported for Cd 
clusters [17]. We note, however, that the experimental ev- 
idence of NMM transitions is rather indirect. Refs. [181 
H9] argue that due to the method of measurement, the 
observed NMM transition may even depend on the for- 
mation process. For a review and discussion, see Ref. [T7] . 

The size range where the NMM transition has been 
proposed is easily accessible by DFT calculations. We will 
contribute to the discussion here by studying the devel- 
opment of HOMO-LUMO gaps as a function of cluster 
size, and by looking at the degree of localization of the 
electronic density. 

Another objective of our work is to prepare for the 
study of metallic clusters in quantum fluid matrices. Tech- 
niques to agglomerate atoms and small molecules in a 
quantum fluid matrix — specifically, in superfluid 4 He — 
have opened a new and versatile way to study the struc- 
tural, electronic and spectroscopic properties of nanopar- 
ticles. The helium droplets can be viewed as ultracold 
nanoscopic reactors, which isolate single molecules, clus- 
ters, or even single reactive encounters at very low tem- 
peratures [20]. Clusters of well-defined composition can 
be formed inside the droplets, and their examination in 
the millikelvin regime has already given important clues 
on magnetism and superconductivity on the nanometer 
scale [2T], 



There is increasing evidence that the presence of a he- 
lium matrix can indeed change the geometric arrangement 
and other properties of the metal cluster ions. For exam- 
ple, cluster growth by capturing in 4 He droplets can lead 
to different isomers, which are not found with other clus- 
ter generation techniques and which may be affected dif- 
ferently by the He environment. Such an effect has been 
observed for clusters bound by hydrogen bridges: Cyclic 
water hexamers were found in helium droplets, but not in 
vacuum [22]. There is also evidence that the feedback of 
the surrounding quantum fluid on the much more strongly 
bound carbon nanotubes is non- negligible; Kim et al. [23] 
point out that this is basically dictated by Newton's third 
law. An unambiguous interpretation of such experiments 
requires understanding of the influence of the 4 He sur- 
rounding the cluster. Also, metallic clusters that are dif- 
ficult to generate in vacuum can be formed in a fluid ma- 
trix [24] . Liquid helium provides an ideal medium for such 
"nanoreactors" because it is transparent in the entire spec- 
tral range from the far IR to vacuum UV, and high spec- 
troscopic resolution, comparable to the gas phase, can be 
achieved. 

Our paper is organized as follows: In Sec. [2] we give a 
brief discussion of the computational methods used. The 
core of our DFT package is a diffusion algorithm for solv- 
ing Schrodinger-like equations; a more extensive analysis 
of the method, including convergence tests, a comparison 
with the implicitly restarted Lanczos method and an as- 
sessment of its performance on parallel computing archi- 
tectures, has been given in Ref. [7]. A separate program 
package for solving the Schrodinger equation for the bound 
states in an arbitrary local potential is also available [25] . 
Sec. [3] turns to the results of our calculations. We briefly 
study Na n and Na+ clusters and determine their structure 
and energetics. A somewhat unexpected feature is that lo- 
cal pseudopotentials predict a distorted ground state con- 
figuration of Na4. We then turn to our discussion of Mg n 
and Mg+ clusters. We present results for the ground state 
configuration, energetics, and stability of these systems 
and examine the appearance and origin of "magic num- 
bers" . The concluding section [4] gives a brief summary of 
our findings. 



2 Computational Methods 

2.1 Spin-Density Functional Theory 

Spin-density functional theory [8] maps the solution of 
the interacting many-electron problem onto that of a non- 
interacting auxiliary system, which is described by the 
Kohn-Sham equations, 



'2m 



V 2 + ^sIK}] (r) 



(r) = £tf (r) , (1) 



^( r ) = EE n il^( r )l 2 = E^( r )- ^ 

a j a 

These are non-linear Schrodinger equations for a set of 
single-particle wave functions ipj (r) in an effective poten- 
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tial Vj£ s [{p a }] (r), where a denotes the spin index and nj 
is the occupation number of the state {j, a}. The effective 
potential 

^ks [{P a }} (r) = ^ext (r, {RJ) + V C [p] (r) + Kc (r) 

(3) 

consists of the external potential ( r ? {Ri}) describing 
the interaction of the valence electrons with the ion cores 
located at the positions {Ri}, the Coulomb term 

V clp}(r) = J ^T-^(rVr\ (4) 

which accounts for the direct electron-electron interac- 
tion, and the exchange-correlation potential Kcdl/^IK 1 *)- 
In this work, the local spin-density approximation [26 is 
used for V^ c , with the correlation functional proposed by 
Perdew and Wang [8]. 

Equations and p]) are solved self-consistently for 
each set of ion positions {R^}. To find the global minimum 
of the total energy, the energy function ^({R^}) must be 
evaluated at many different points in configuration space. 
Thus, an indispensable prerequisite for performing an un- 
restricted minimization is a fast and reliable method to 
solve the Kohn-Sham equations. 



2.2 Diffusion algorithm 

In recent years, real-space methods for solving the Kohn- 
Sham equations 0, ([2| have become increasingly popular 
[27tl28]. They are easy to implement, free from modeling 
uncertainties, and well suited for modern massively par- 
allel computers. The calculations in this work have been 
performed using the program package limerec [4], which 
solves the Kohn-Sham equations on a real- space grid using 
an efficient diffusion algorithm: The lowest n eigenstates 
of the one-body Schrodinger equation 

{T + V)^{v) = E^{v) (5) 

are obtained by repeatedly applying the imaginary time 
evolution operator 

T ( e ) = e -<T+V) (6) 

to a set of trial functions {^(r)}, which are orthogo- 
nalized after each step. Above, T is the one-body kinetic 
energy operator, and V the potential. The operator T(e) 
multiplies states corresponding to higher energies by expo- 
nentially decreasing weights. It thus gradually filters out 
high-energy states, making the procedure converge to the 
lowest n eigenstates of the operator T + V. 

The evolution operator (|6| can not be calculated ex- 
actly. A simple approximation is the so-called split oper- 
ator form [29] . 

r2 ( e ) = e -htV e -eT e -yv = e -e(T + V + 0(S)) . (?) 

it is accurate up to second order in the timestep e. 



Factorizing T(e) obviously reduces the problem of cal- 
culating the exponential of the full Hamiltonian to the 
problem of dealing with its parts separately. For a local 
operator V^(r), the action of e~2 ey ( r ) on a function ^-(r) 
is just a vector- vector multiplication. The operator e~ eT 
is local in momentum space, its action on a given state 
needs one Fast Fourier Transform (FFT) pair (forward 
and backward) in addition to a vector-vector multiplica- 
tion in Fourier space. If ^j( r ) is represented by its values 
at N grid points, these operations scale with O(N) and 
0(7Vlog7V), respectively. 

The timestep-dependent error 0(e 2 ) introduced by the 
decomposition process imposes an upper limit for the time 
step if a certain accuracy has to be reached. On the other 
hand, the whole procedure needs fewer iterations when 
the time step is large. To achieve faster convergence, it 
would be desirable to use higher order algorithms at larger 
time steps. Unfortunately, Sheng [30] and Suzuki [31 have 
shown that no factorization of the evolution operator ^ 
into a product of the operators e~ eaV and e~ e ^ T beyond 
second order can have all positive coefficients a,/3. This 
so-called forward time step requirement is, however, essen- 
tial for the numerical stability of the algorithm. Fourth 
order algorithms with exclusively positive time steps have 
been derived by Suzuki [32,33 and Chin [34 by intro- 
ducing an additional correction to the potential of the 
form [V, [T, V]]. This correction term has first been used by 
Takahashi and Imada [35] . We have shown previously [36] 
that these forward fourth-order algorithms can achieve 
similar accuracy at more than an order of magnitude larger 
step sizes compared to second-order splitting (P7|). 



2.3 Multi-product expansion 

Recently, we have made important progress in deriving 
higher-order factorization schemes [7,25] by using a linear 
combination of second-order propagation steps, 

n 

T 2 „ (e) = J2 (I) = e-(™^ 2 "» . (8) 

k=l 

The coefficients Ck are given in closed form for any k 
[37] . We have recently implemented a diffusion algorithm 
based on such an expansion [7|[25] . One 2n-th order prop- 
agation step 72n( e ) requires n(n+l)/2 second-order prop- 
agation steps. This is more than outweighed by the fact 
that much larger timesteps can be used for the higher- 
order algorithms due to the smaller time-step error. For 
instance, we have shown in Refs. [7)125] that the 12-th 
order algorithm, which requires 21 second-order propaga- 
tions, yields the same accuracy as the underlying second- 
order scheme at a time step about a factor of 1000 larger. 
The advantage of high-order propagation methods is par- 
ticularly compelling on parallel computer architectures: 
The calculation of the action of the evolution operator ([6| 
on the wave functions can be parallelized efficiently by 
simply distributing the wave functions ^j( r ) across dif- 
ferent processors. As higher-order algorithms need much 
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fewer, but more expensive propagation steps, the relative 
weight of orthogonalization, which can be parallelized less 
efficiently, decreases. Consequently, the parallel speedup 
ratio is higher for higher-order algorithms [7]. 

2.4 Pseudopotentials 

Chemical binding properties are almost exclusively deter- 
mined by the valence electrons, a fact that is particularly 
true for metals like sodium and magnesium where the s- 
type outer electrons are well separated from a noble gas 
type ion core. The idea to ignore the strongly bound core 
electrons in calculations and to reduce the nucleus and 
the core electrons to a "black box", an ion core that in- 
teracts with the valence electrons by means of an effective 
pseudopotential, dates back to Fermi [38] . 

Pseudopotentials used in modern electronic structure 
calculations generally fall into two categories: One type, 
empirical pseudopotentials, uses an analytic model po- 
tential with parameters that are fitted to experimental 
data. So-called ab-initio pseudopotentials, on the other 
hand, are obtained by inverting the free-atom Schrodinger 
equation for a given reference configuration [39], and en- 
force the pseudo-wave functions to coincide with the true 
all-electron wave functions outside a given core radius. 
The pseudopotentials in this group are usually non-local. 
The diffusion algorithm, which was presented for local po- 
tentials in Sec. |2.2[ has also been implemented for non- 
local pseudopotentials [6 . The program package limerec 
supports both local as well as norm-conserving non-local 
pseudopotentials in the form generated by the fhi98pp 
[40] pseudopotential generator. Most calculations in this 
work have been carried out using the local, empirical po- 
tentials of Fiolhais et al. [9 . For the purpose of testing the 
sensitivity of our results to the pseudopotential used, we 
have performed the calculations for Na clusters also with 
the potential proposed by Kummel et al. [10]. In one case, 
the local potentials predicted a ground state structure dif- 
ferent from what has been reported in the literature. In 
this case, we have checked the calculations using non-local 
Troullier-Martins pseudopotentials [TT] , 

2.5 Structure optimization 

For structure optimization, we have used a simulated an- 
nealing procedure based on Langevin [4Tj dynamics of the 
ions and a steepest descent method. 

The Langevin method is similar to the Metropolis al- 
gorithm [42] , but additionally takes into account the influ- 
ence of forces on the particles. The ions at positions {Ri} 
are moved according to 

R (^) = R W + fag(i) _^ F W, (9) 

where g(i) is a sequence of Gaussian random numbers, Sx 
is a factor that determines the amplitude of the random 
walk of the atoms, is the force on the i-th atom in 



simulation step N, and j3 = l/(fc j eT), where T is an artifi- 
cial temperature. The forces on the ions can be calculated 
from a single DFT calculation at fixed ion positions {R^} 
using the Hellman-Feynman theorem [43 ] l44 ] . After each 
move, the energy difference 

AE = E[{R ( i N+1) }]-E[{-R! l N) }} (10) 

is calculated, where £?[{R^}] is the energy of the clus- 
ter at fixed ion positions {Hi} at step N, as calculated 
by DFT. If AE < 0, the move is accepted in any case, 
for AE > 0, the move is accepted with a probability 
exp(— f3 AE). Compared to the Metropolis algorithm, Lan- 
gevin dynamics allows for higher acceptance rates. 

Langevin simulated annealing is guaranteed to con- 
verge to the global energy minimum asymptotically, i.e., 
in the limit of infinitely many annealing steps. For a de- 
tailed analysis of the convergence properties of the meth- 
od, including the finite-time behavior, see Ref. [45]. In 
practical calculations, cooling the system faster or termi- 
nating the annealing process earlier leads to an increased 
probablility of freezing it in a higher-lying, local energy 
minimum. One way to verify the robustnes of the proce- 
dure is to repeat it from different starting configurations. 
To find the ground-state configurations and isomers re- 
ported in this work, we have thus proceeded as follows: 

1. We have generated a number of random starting con- 
figurations of ion core positions for each cluster, only 
ensuring that the distances between the ion cores are 
large enough to avoid divergences. 

2. For each starting configuration, we have performed a 
simulated annealing using Langevin molecular dynam- 
ics: After a suitable thermalisation period, the temper- 
ature in Eq. ([9| was slowly decreased to "freeze" the 
cluster in a low- lying minimum of the energy surface. 

3. As soon as the annealing procedure had settled close to 
a local minimum at low temperature, we switched to 
a steepest descent method to accurately determine the 
position of the minimum. This was implemented by 
using the Verlet molecular dynamics method for the 
ion cores and setting, after each step, all ion velocities 
to zero. This changed the ion positions by typically 
about 0.1 clq. 

4. If only ground- state configurations are of interest, the 
annealing process can in principle be stopped as soon 
as it becomes clear that it converges towards a higher 
energy than a previously found configuration. To find 
isomers, we have nevertheless completed the full an- 
nealing procedure for all Na n and Na+ clusters shown 
here, as well as for Mg n and Mg+ up to n = 15. 

5. All found configurations were checked for congruency, 
i.e., for being identical up to rotations and/or transla- 
tions. It turned out that all structures with an energy 
difference greater than 10 -4 Ry per electron were iso- 
mers. In four cases we have also found incongruent 
configurations with an energy difference of less than 
10 -4 Ry per electron. Two of these will be discussed 
below. 

To estimate the discretization error, we have also cal- 
culated, for some sample clusters, the total energy while 
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rigidly moving and rotating the ions in different directions 
in the simulation box. For all clusters, the discretization 
error is around 10 -5 Ry, which is certainly much smaller 
than the uncertainty implicit to the local density approx- 
imation. 

In order to find a locally stable configuration, we have 
typically carried out about 1000 Langevin moves, slowly 
decreasing the artificial temperature T = 1/(&b/3). These 
were followed by several hundred steepest descent moves. 
Therefore, each locally stable configuration shown here in- 
volved between 1000 and 2000 LDA calculations for the 
cluster under consideration. The diffusion algorithm de- 
scribed above is particularly efficient for such a strat- 
egy, because it iteratively improves a given set of initial 
wave functions, and thus profits from good starting val- 
ues. Since the ions are moved only very little during one 
annealing step, the wave functions do not change much, 
and the diffusion algorithm for subsequent steps converges 
within just a few iterations. 
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Fig. 1. (color online) The binding energy per atom of neu- 
tral Na n clusters compared with earlier work as cited in the 
legend. The results of this work have been obtained using the 
potentials of Fiolhais et al. [9] and those by Kiimmel et al. [10] . 



3 Results 

3.1 Na and Na + clusters 

Na n and Na+ clusters are the most frequently studied type 
of metallic clusters, for a recent review see Ref. [T2] . We 
have examined these systems basically to highlight the 
differences between Na n and Mg n clusters. We have per- 
formed (S)DFT calculations for small Na n and Na+ clus- 
ters up to n = 10. The calculations were carried out on 
a regular, Cartesian real-space grid with 64 3 mesh points 
and a resolution of 0.5 ao, where ao is the Bohr radius. 
We have determined that this is the coarsest discretization 
that is permitted to get reliable results. Our discretization 
corresponds to an energy cutoff of 537 eV in plane- wave 
calculations. 

All calculations in this section have been carried out 
using the empirical pseudopotentials of both Fiolhais et 
al. [9] and Kiimmel et al. [10 , for the purpose of testing 
the sensitivity of the results to the pseudopotential used. 
In one case, we have additionally used norm-conserving 
non-local pseudopotentials of the Troullier-Martins type 
[U, see Fig. [4] 

3.1.1 Energetics 

The binding energy for neutral and singly charged clusters 
is defined by 

where E\ and E n are the ground state energies of a single 
Na atom and a Na n cluster, respectively, and E^ and E+ 
are their counterparts for the charged clusters. 

Figs. [I] and [2] show the binding energy per atom for 
Na n and Na+ clusters in comparison with previous work 




Cluster size n 

Fig. 2. (color online) Binding energy per atom of singly ion- 
ized Na+ clusters compared with earlier work as cited in the 
legend. The results of this work have been obtained using the 
potentials of Fiolhais et al. [9] and those by Kiimmel et al. [10] . 



pT2 |l46|l4T| l48] . where we find generally good agreement. 
For neutral clusters we have found that the pseudopoten- 
tial proposed by Fiolhais et al. [9] and that of Kiimmel et 
al. [10] give almost identical results. For charged clusters, 
the difference in the binding energy due to different pseu- 
dopotentials is comparable to the discrepancy between dif- 
ferent calculations. 

Consistent with earlier work, we observe a particularly 
high stability for Na2, Nag, NaJ and Na^, where the bind- 
ing energy per atom is higher than the one of the cluster 
with one additional atom. This is known to be due to 
electronic shell closure and predicted by ellipsoidal jellium 
models [49]. We also observe an additional peak for Na^~ 
which can not be explained by shell closures. This peak 
is consistent with earlier theoretical calculations [50]; it 
has also been observed in Na+ fragmentation experiments 
[51]. It is the onset of an even-odd alternation which has 
been attributed to electron pairing effects [52] . 
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Fig. 3. (color online) Ground state configurations of neutral 
Na n and singly charged Na+ clusters up to n — 10. The results 
reported here were obtained using the potential of Fiolhais 
et al. 9 1. Note that we show, in the case Na4, the rhomboid 
isomer. Local pseudopotentials predict an additional distorted 
configuration slightly lower in energy, see Fig. [4] Eb, n /n is the 
binding energy per atom, (d) n the average nearest neighbor 
distance, and d m in,n and d max ,n are the smallest and largest 
nearest neighbor distances, respectively. 



3.1.2 Geometric structure 

An overview of the annealed ground state structures for 
Na n and Na+ is given in Fig. pi 

The agreement with earlier work is consistent with 
what is expected from the energetics. The only noticeable 
deviation from literature has been observed for Na4, where 
the energy minimization with both local pseudopotentials 
converges towards a distorted rhomboidal structure, see 
Fig.[4j The symmetric rhomboidal structure, which is usu- 
ally reported as the ground state in the literature (see, for 
instance, Ref. [48 ) appears as an isomer only 3.96- 10 -4 Ry 
higher in energy. When using non-local pseudopotentials 
of the Troullier-Martins type [11] , the distorted structure 
ceases to be a local minimum, and the symmetric configu- 
ration is found for the ground state. The effect thus seems 
to be connected to the locality of the potential. 



<$> <I> 

Fig. 4. (color online) Ground state configurations of Na4. Left 
pane: "distorted" rhomboidal configuration which is found as 
the ground state with both local pseudopotentials. Right pane: 
Symmetric configuration, found as slightly higher-lying iso- 
mer with local potentials (energy difference: +3.96 • 1(T 4 Ry). 
For non-local Troullier-Martins pseudopotentials, the distorted 
configuration is not stable, and the symmetric configuration 
shown right is the ground state structure. 



3.1.3 Kohn-Sham energies 

The Kohn-Sham energies of the ground state structures 
for Na n and Na+ are depicted in Figs. [5] and [6j For odd 
electron numbers we have used the LSDA, which implies 
that electrons with different spins can have different wave 
functions and energy eigenstates. The lowest states in these 
figures correspond to the ls-states. The lp-states are sep- 
arated by a pronounced gap and can be clearly identified, 
as well as the onset of the ld-states. 

The energy levels show a good overall agreement with 
the single-particle levels predicted by the Clemenger-Niels- 
son model [53] . see Fig. 5(a) in Ref. [49] . This indicates 
that ellipsoidal jellium models give a good qualitative de- 
scription of Na clusters. One noteable difference is Nag: 
While clusters with closed shells usually are spherically 
symmetric and thus show degenerate p-states in jellium 
models, this degeneracy has been lifted in our calculation. 
This is a result of the symmetry breaking due to the ion 
cores, and has been observed previously, see, e.g., Fig. 5(d) 
in Ref. [49]. 

The Clemenger-Nilsson model predicts shell closures 
for Na2, Nag, Na^ and Nag" [49,53 . Correspondingly, the 
neighboring clusters show a significantly decreased gap in 
at least one spin channel in Figs. [5] and [6] In addition to 
these shell closures, we also observe the even-odd alterna- 
tion attributed to the electron pairing effects mentioned 
above [52] . Due to this effect, also Na6 and Na^" show a 
particularly large gap. 



3.2 Mg and Mg + clusters 

The focus of our work is the calculation of structures and 
energetics of Mg n and Mg+ clusters. We have again used 
the procedure outlined in Sec. 2.5 to find ground states 
and low-lying isomers of Mg n and Mg+ clusters up to 
n = 30. The ion cores of the Mg atoms have been described 
by the local potentials of Ref. [9] . 

Due to the shorter range of the Mg pseudopotentials 
compared to those for Na, all calculations have been per- 
formed on a regular, Cartesian real-space grid with 96 3 
mesh points and a resolution of 0.375 ao- This corresponds 
to an energy cutoff of 954.5 eV in a plane- wave calculation. 
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Fig. 5. (color online) Kohn-Sham energy levels of the neutral 
Na n clusters as a function of cluster size. Left pane: Energy 
levels for spin- up electrons. Right pane: Energy levels for spin- 
down electrons. The HOMO levels are marked by the lower 
(blue) solid line, the LUMO levels by the upper (red) solid line. 
The even electron numbers have been calculated with LDA, the 
odd ones with LSDA, hence the energy levels of spin-up and 
spin-down electrons of the configurations with even electron 
numbers are the same. 



U P down 
-0-1 | 1 1 1 1 1 1| 1 1 1 — 




2468 10 2468 10 

Cluster size n 

Fig. 6. (color online) Kohn-Sham energy levels of the singly 
ionized Na+ clusters as a function of cluster size. See Fig. H] for 
further explanations. 



A comment is in order concerning the expected accu- 
racy of the applied theory. This can be assessed by com- 
paring predictions for small systems with results obtained 
by more elaborate methods. Refs. [18] and [54], for in- 
stance, compare the bond lengths of Mg3 and Mg4 ob- 
tained from coupled cluster (CC) and DFT calculations. 
The CC calculations included single and double excita- 
tions as well as perturbative triples [CCSD(T)], the DFT 
calculations were all-electron calculations using the BP86 
correlation functional. The results of both methods agree 
with ours within a few percent. This is certainly better 
than expected because the pseudopotentials used in our 
work have been designed to reproduce bulk properties. 
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Fig. 7. (color online) Binding energy per atom, Eb, n /n, of Mg n 
clusters in comparison with earlier calculations as cited in the 
legend. 



3.2.1 Energetics 



To safely bracket the regime where a non-metal to metal 
(NMM) transition has been proposed [14,15 , we have cal- 
culated structures up to n = 30. The binding energy for 
neutral and singly charged clusters is defined in Eq. ( pT| , 
where E\ and E n are the ground state energies of a single 
Mg atom and a Mg n cluster, and and E+ their coun- 
terparts for the charged clusters, respectively. By repeat- 
ing the annealing procedure with different initial configu- 
rations, we have also calculated isomers for ion numbers 
n < 15; for larger clusters the number of isomers grows 
rapidly and identifying all of them becomes extremely 
time consuming. 

Figures [7] and [8] show our results for the binding en- 
ergy per atom of the Mg n and Mg+ clusters, compared to 
a number of previous calculations [18,54,55,56,57 . The 
deviations between the calculations are due to the use 
of different energy functionals and pseudopotentials, and 
give an idea of the expected accuracy of the underlying 
physical model. The most apparent discrepancy between 
our results and earlier work is the greater smoothness of 
the binding energy as a function of the number of atoms. 
This is likely due to our elaborate annealing procedure in 
combination with high spatial resolution. 

For a quantitative assessment of cluster stability, we 
turn our attention to the possible fragmentation channels. 
A Mg n cluster can dissociate into two smaller clusters of 
sizes n — m and m. The fragmentation energies are defined 
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Fig. 8. (color online) Binding energy per atom, E^ n /n, of 
Mg+ clusters in comparison with earlier calculations as cited 
in the legend. 



AE n = max [E h>n _ m 

l<m<n ' 

AE+= max [E+ n _ m 

1 < rn < n ' 



(12) 



If this quantity is negative, the cluster is stable against 
spontaneous fragmentation. Nevertheless, for each clus- 
ter, we denote the channel with the highest fragmentation 
energy as the preferred fragmentation channel, even if the 
cluster is stable. Figs. [9] and [To] show the dependence of 
the fragmentation energy on cluster size. We found that 
all Mg and Mg + clusters are stable. The preferred frag- 
mentation channel for neutral clusters is always 



Mg n — ► Mg„_ 1 + M gl . 
For singly ionized clusters it is, correspondingly, 
Mg+ — ► Mg+_! + M gl 



(13) 



(14) 



for all cluster sizes. This is due to the monotonic increase 
of the binding energy with the number of atoms, in con- 
trast to Na clusters, which often decay as Na n — » Na n _ 2 + 
Na 2 . 

Another criterion for cluster stability [49] involves the 
second energy difference, 



A>M+> EE E 



(+) 



+ E, 



(+) 

n+l ' 



(15) 



where the superscript (+) indicates that the formula ap- 
plies to both neutral and charged clusters. Eq. ( 15 ) can be 
understood as a discretized second derivative oTthe elec- 
tronic energy with respect to the ion number. If A 2 E < ^ ) > 

0, then Mgn + ^ is stable with respect to disproportionation 
of two Mgn + ^ clusters into Mg^} 1 and Mg^_\. A higher 
value indicates a more stable cluster. The evolution of 
these quantities with cluster size among some represen- 
tative results found in the literature are shown in Figs. 11 
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Fig. 9. (color online) Fragmentation energies of Mg n clusters 
in the energetically preferred channel, AE n , for neutral Mg n 
clusters up to n = 30, see Eq. (12). Results of some earlier 



works are shown as cited in the legend. 
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Fig. 10. (color online) Fragmentation energies of Mg^ clusters 
in the energetically preferred channel for singly ionized Mg^ 
clusters, AE^ , up to n — 30, see Eq. (12). Results of some 



and 12 , respectively. 



earlier works are shown as cited in the legend. 



We can use both the fragmentation energies and the 
second energy differences to identify extraordinarily sta- 
ble clusters. To find those, we look for clusters that show 
a distinct minimum in AE^ and a distinct maximum in 
in the region where A2E^ > 0. This procedure 
suggests that neutral Mg n clusters with 4, 7, 10, 17, 21, 23 
and 27 atoms are particularly stable. For the ionized Mg+ 
clusters, the most stable sizes are 5, 8, 10, 14, 20, 23 and 
27, respectively. Diederich et al. [15] have measured the 
abundance of Mg+ clusters grown in ultracold liquid he- 
lium nanodroplets using mass spectroscopy. They report, 
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Fig. 11. (color online) Second energy differences of neutral 
Mg n clusters, A^E n ^ up to n — 29, see Eq. (15). A number of 



different results reported in the literature is shown. 
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Fig. 12. (color online) Second energy differences of singly ion- 
ized Mg+ clusters, A^E^ , up to n — 29, see Eq. ( 15 ). A number 



of different results reported in the literature is shown. 



in the area of cluster sizes covered here, distinct maxima 
in the abundance distribution at clusters sizes 5, 10 and 
20, as well as a local maximum at 8. All of these "magic 
numbers" are consistent with the stable cluster sizes pre- 
dicted by the above procedure. 

The clusters Mg^" 4 , Mg^ and Mg^, which are also pre- 
dicted to be particularly stable by our analysis, do not 
show up as abundancy maxima in the experiment. Two of 
these clusters, Mg+ 4 and Mg^, seem to have a somewhat 



peculiar geometrical structure, see Sec. 3.2.2 



The evolution of the adiabatic ionization energies, 

^ion,n = ^ - E n , (16) 
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Fig. 13. (color online) Adiabatic ionization energies, i?i on ,n, 
of Mg clusters up to n = 30 in comparison with earlier calcu- 
lations as cited in the legend. 



with cluster size is shown in Fig. 13 The almost mono- 
tonically decreasing values are in good agreement with 
previous calculations [55| f56] . 



3.2.2 Geometric structure 

Fig. [14] shows the ground state geometries of the neutral 
and charged Mg n clusters up to n = 11. Displaying the 
individual ion positions for larger clusters does not provide 
much useful information, these positions can be found in 
Ref. [5H]. 

To examine the structure of larger clusters, we look 
at the pair distribution function of the ions. The classical 
two-body density of the ion cores, 



p 2 (r,r')^^«5(r 



R)<5(r'-R,), 



(17) 



contains the most useful information about the geometric 
structure of the cluster. To reduce the number of vari- 
ables, we rewrite P2(r, r') in terms of the relative coordi- 
nate x = r — y' and the center of mass coordinate and 
average over the center of mass coordinate as well as over 
the remaining angular degrees of freedom. This defines the 
so-called average radial distribution function g n (x), 



9n{x) 



1 



47rx 2 pon 



IR-R.-I), (18) 



where po is the average bulk density. This procedure yields 
the exact pair distribution function for homogeneous sys- 
tems. Finally, to obtain a smooth function, we have broad- 
ened the distribution of pair distances x by a Gaussian of 
width e, 



9nA X ) 



47ry27rer 2 pon 



^3 



.Or-lRi-R,-!) 2 



(19) 



For small e, Eq. ([19]) is identical to the discretized expres- 
sion in Eq. (fl8l). 
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Fig. 14. (color online) Ground state configurations of neu- 
tral Mg n and singly charged Mg+ clusters up to n = 11. The 
ground state structures of all Mg n clusters with n < 30 can 
be found in Ref. 



[58] ; isomers will be discussed in Sec. 3.2.5 
Eb,n/n is the binding energy per atom, (d) n the average near- 
est neighbor distance, and d m in,n and d ma x,n are the smallest 
and largest nearest neighbor distances, respectively. 



The peaks in g n ,e(%) correspond to a high probability 
of finding two atoms at a distance x and, in the bulk 
limit, correspond to the sequence of n-th nearest neighbor 
distances £ n . For face-centered cubic (fee) and hexagonal 
close packed (hep) structures, the ratios between first and 
second, and first and third nearest neighbor distances are 

6 : & = V2, £3 : £x = Vs (fee and hep). (20) 

The pair distribution functions g n ^ e (x) for a few selected 
neutral and singly ionized Mg clusters are shown in Figs. 15 
and 16 While Mg 30 is still far away from the bulk limit, 



one can nevertheless clearly identify a sequence of peaks 
in the plot of g3o je (x). I n this case, we have determined the 
positions of these peaks to be £1 = 6.275 ao, £2 = 8.375 ao 
and £3 = 10.725 ao, respectively. The corresponding ratios 
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Fig. 15. (color online) Pair distribution function for selected 
Mg n clusters. The inset shows the distribution functions for 
x > 7.5 ao magnified (right scale). We have chosen a broaden- 
ing parameter of e = 0.15 ao- 
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Fig. 16. (color online) Pair distribution function for selected 
Mg^ clusters. The inset shows the distribution functions for 
x > 7.5 ao magnified (right scale). We have chosen a broaden- 
ing parameter of e = 0.15 ao. 



of peak positions are 
6:6 = 1-34 (-5.24%), 



6:6 = 1-71 (-1.27%), (21) 



where the values in parentheses give the error relative to 
the known ratios for the bulk fee and hep lattices as given 
in Eq. ( |20| ) . From our calculations one would therefore con- 
clude that magnesium clusters approach, with increasing 
cluster size, either fee or hep structure. The Mg 30 cluster 
is still too small to reliably extract a fourth peak from 
^30,e(^), which would allow to distinguish between these 
two structures. However, up to the cluster sizes taken into 
account here, our results are consistent with the known 
hep symmetry of bulk magnesium. 

We finally draw attention to another remarkable piece 



of information that can be gained from Figs. 15 and 16 
When comparing the average pair distribution functions 
for the different clusters in the range from 8ao to 10 ao, 
it is evident that g n , e (x) is zero for the clusters Mg^" 4 , 
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Fig. 17. (color online) Average nearest neighbor distance (d) n 
for neutral Mg n and charged Mg+ clusters as a function of 
cluster size n. Also shown is the bulk limit (dotted line) and 
the theoretical results of Refs. [55] and |56| . 



Mg23 and Mg^, while it is non-zero for all other clusters 
examined here. It can be concluded that those clusters do 
have a geometric structure that is qualitatively different 
from that of the other "magic" clusters. 

From the geometric configurations of the Mg n clus- 
ters we can calculate the average nearest neighbor dis- 
tance (d) n and the smallest and largest nearest neighbor 



distances d n 



and d rr 



respectively. The numerical 



values for these quantities up to a cluster size of n = 11 
are given in Fig. 14 The evolution of (d) n with cluster 



size for both neutral and singly ionized clusters is shown 
in Fig. 17, together with results from Refs. [55 and [56] . 
The dotted line in this figure marks the bulk limit of 
(d) = 6.07 ao for solid magnesium. Evidently, at n — 30 
both neutral and charged clusters are already close to that 
limit. 



3.2.3 Kohn-Sham energies 

The Kohn-Sham levels of the ground state structures of 
and Mg+ are shown in Figs. [l8| and 



M 



19 



For the Mg+ 



clusters, the LSDA has been used; we therefore show the 
energy levels for both spin channels in separate plots. 

In contrast to the HOMO-LUMO gaps of the sodium 
clusters shown in Figs. [5] and [6| the gaps of the magnesium 
clusters do not show strong signs of electronic shell clo- 
sures: The gap decreases rather monotonically with cluster 
size; the pronounced "zigzag" structure of the HOMO and 
LUMO levels found in sodium is missing. This is an indi- 
cation that small Mg clusters are more covalently bound 
and electrons are more localized in comparison to sodium 
clusters. This prediction is confirmed by the respective 



electron densities, see Sec. 3.2.4 



In a bulk system, the closing of the HOMO-LUMO 
gap can be considered as a sign of an NMM transition. 
In finite-sized systems like clusters, the gap between en- 
ergy levels is always finite, and the answer to the ques- 
tion whether a material is a metal or not depends on 
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Fig. 18. (color online) Kohn-Sham energy levels of the neutral 
Mg n clusters as a function of cluster size. The HOMO levels 
are marked by the lower (blue) solid line, the LUMO levels by 
the upper (red) solid line. 



2 



o 



o 



3 

o 




15 20 

Cluster size n 



Fig. 19. (color online) Kohn-Sham energy levels of the 
ionized Mg^ clusters as a function of cluster size. Top 
Energy levels for spin-up electrons. Bottom pane: Energy 
for spin-down electrons. 
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Fig. 20. (color online) Contour plot of the electron density of Fig. 21. (color online) Contour plot of the electron density of 



Na7 cluster in the symmetry plane, in units of 10 3 ao 3 . 



a Mg7 cluster in the symmetry plane, in units of 10 3 ao 3 . 



its temperature [17] . Specifically, a temperature of 300 
Kelvin corresponds to roug hly 1.9-10- 3 Ry. The calculated 



HOMO-LUMO gaps shown in Figs. [18] and [19] are on the 
order of 0.02-0.10 Ry. Therefore, Mg clusters with up to 
30 atoms are not metallic by this criterion. Note that the 
LDA is known to systematically underestimate band gaps 
(see, e.g., Ref. [59 J, and is thus erring on the side of cau- 
tion when predicting non-metallicity. 

At first glance, this conflicts with Refs. [HJ 15, 16 , where 
the NMM transition has been observed for Mg clusters in 
helium at around n = 20. However, Refs. [T8]H9] suggest 
that the observed NMM transition most likely refers to 
negatively charged Mg~ clusters used in these experiments 
and thus does not disagree with our findings. 




3.2.4 Electron density 

The most telling documentation of the difference between 



Na n and Mg n clusters is the electron density. In Figs. 20 



and 21 we show, as examples, contour plots of the electron 
density of the Na7 and Mg7 ground state configuration 
in the symmetry plane. We have chosen these examples 
because their ionic configuration is very similar, see Figs. [3] 
and El 



From Figs. 20 and 21 it is clear that the electrons in 
a Na7 cluster are much more delocalized than in a Mg7 
cluster. If there were an NMM transition in the range of 
cluster sizes under consideration here, one would expect 
a derealization of electrons and a much broader electron 
density for larger clusters. In Fig. [22] we show therefore 
a contour plot of the electron density in a Mg3o cluster. 
Note that the scales are the same as in Fig. |2l] Evidently, 
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Fig. 22. (color online) Contour plot of the electron density of 
a Mg3o cluster in the 2/2-plane, in units of 10 -3 ao -3 . 



the electrons are just as strongly localized as in Mg7, indi- 
cating that neutral clusters have, at that ion number, not 
yet undergone an NMM transition. 



3.2.5 Isomers 



Tables 1 and 



i 



show the binding energies per atom, 



(+) 



for the grouna states of neutral Mg n and singly charged 
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Mg n 




6[10- 3 Ry] 


c[10- 3 Ry] 


d[10- 3 Ry] 


1V± &2 


004904 








1V± &3 


0.011110 


-3.653 






iV1 &4 


020072 








1V± &5 


021798 








iV1 &6 


023444 


-0.291 






Mg 7 


0.026227 


-0.450 


-0.606 




Mg 8 


0.028047 


-0.278 


-0.798 


-1.215 


Mg 9 


0.031079 


-1.360 






Mg 10 


0.033874 


-0.357 


-1.827 




Mg n 


0.034759 


-1.572 


-2.004 




Mg 12 


0.034859 


-0.547 


-0.762 




Mg 13 


0.035357 


-0.226 


-0.447 


-1.219 


Mg 14 


0.037420 


-0.563 


-0.590 


-0.889 


Mg 15 


0.039016 


-0.869 







Table 1. Binding energies per atom, Eh, n /n, of the ground 
state configurations and a few isomers of neutral Mg n clus- 
ters up to n = 15. The second column shows Eh, n /n for the 
ground-state configuration, columns (b)-(d) show the differ- 
ence in binding energy per atom between the ground state and 
the second, third, and fourth isomer. 



Mg+ 




6[10- 3 Ry] 


c[10- 3 Ry] 


d[10- 3 Ry] 


Mg+ 


0.048724 








Mg 3 + 


0.053346 








Mg+ 


0.052536 


-1.382 






Mg+ 


0.053102 


-0.902 






Mg 6 + 


0.052759 


-1.296 


-1.622 


+2.235 


Mg+ 


0.051013 


-0.013 


-0.182 




Mg 8 + 


0.051870 


-0.008 


-0.892 




Mg+ 


0.051520 


-0.024 


-0.214 




Mg+ 


0.053659 


-0.049 


-0.485 




Mg+ 


0.053859 


-0.333 


-1.496 




Mg^ 2 


0.053125 


-0.304 


-1.051 




Mgi 


0.053817 


-0.184 


-0.379 


-0.425 


Mg+ 


0.054982 


-0.335 


-0.375 


-0.577 


Mg+ 


0.055399 


-0.268 







Table 2. Binding energies per atom, E^ n /n, of the ground 
state configurations and a few isomers of charged Mg+ clus- 
ters up to n — 15. The second column shows E^ n /n for the 
ground-state configuration, columns (b)-(d) show the differ- 
ence in binding energy per atom between the ground state and 
the second, third, and fourth isomer. 



Mg+ clusters, together with the differences in binding en- 
ergy per atom of identified isomers to the corresponding 
ground state. 

A remarkable feature of these results can be seen in Ta- 
ble^ For all clusters between Mg^ and Mg+ , the lowest- 
lying isomer is energetically almost degenerate with the 



ground state. Figs. [23] and [24] show the corresponding ge- 
ometric structures for Mg^ and Mgg", respectively. In the 
case of Mg^ , see Fig. 



23 



both the ground-state and the 
lowest-lying isomer are composed of the same Mgg" core 
structure, with the additional atom added at two differ- 
ent positions. For Mgg", shown in Fig. 24 , the situation 



Fig. 23. (color online) The two structures of Mg+ with 
lowest energy. Left pane: ground state configuration. Right 
pane: lowest-lying isomer (energy difference: —1.3 • 10 -5 Ry, 
see Tab. [5). 





Fig. 24. (color online) The two structures of Mgg with 
lowest energy. Left pane: ground state configuration. Right 
pane: lowest-lying isomer (energy difference: —2.4 • 10 -5 Ry, 
see Tab. [2]). 



seems to be more complicated. Here the atom positions 
and bonding angles for the two configurations are consid- 
erably different. 



4 Conclusion 

We have presented a systematic study of neutral and singly 
charged Na n and Mg n clusters. We have employed a real- 
space algorithm that is particularly well suited for the 
problem at hand: It allows a fine coordinate-space res- 
olution that is able to deal with even small annealing 
steps. The process is iterative and profits from the knowl- 
edge of accurate solutions of near-by configurations. We 
found, during the annealing process, that the HOMO- 
LUMO gaps are indeed very sensitive to the details of 
the ground state structure. Thus, the last annealing steps 
were of the order of 10 -4 ao- We have verified the numeri- 
cal stability of our procedure by moving and rotating the 
cluster with respect to the coordinate space mesh. The to- 
tal energy had uncertainties comparable to what is gained 
during the final annealing steps, but the relative energies 
between neighboring configurations were independent of 
the precise location and orientation of the cluster. 

The relatively strongly localized electron density in Mg 
clusters also suggests that a simple jellium model can not 
describe the cluster properties accurately. In contrast to 
alkali clusters in general, where the delocalized electrons 
dominate the cluster structure, the geometric symmetry 
seems to be more important for Mg clusters. The highly 
symmetric Mgn cluster has, for instance, only a slightly 
higher binding energy per atom than Mgio, for which the 
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jellium model predicts an electronic shell closure. Never- 
theless, the stability analysis depicted in Fig.[9]shows that 
Mgio is still more stable than Mgn. This is different for 
Mg2o and Mg2i: The jellium model predicts a shell closure 
for Mg2o, but our calculations show that Mg2i is slightly 
more stable than Mg 2 o- 

Because Mg n and Mg+ clusters with up to 30 do not 
show closing HOMO-LUMO gaps, we conclude that these 
clusters are not metallic. This is consistent with the find- 
ings of Refs. [18 , 19] that the observed NMM transition for 
Mg clusters in helium at around n = 20 [14,15 ,16 most 
likely refers to the negatively charged Mg~ clusters used 
in the corresponding experiments. 
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